C.................... RSTEMC_EXB.FOR .............................
C.... this subroutine sets up the error functions of the time-dependent
C.... ion and electron energy equations for subsequent solution by the 
C.... Newton iterative procedure. It calls a function for Ti and one for
C.... Te. This version was created by Silvy John to implement the
C.... flux preserving method (1996). The flux preserving approach
C.... is explained in the paper by Torr et al., JGU, December 1990
      SUBROUTINE TFIJ(J,ILJ,IPR,N,TI,F,JB,JBS,V,HE)
      IMPLICIT DOUBLE PRECISION(A-H,K-L,N-Z)
      DIMENSION N(4,401),TI(3,401),F(20),V(2),HE(401)
      COMMON/VN/U(2,401),BG(401),BM(401),GR(2,401),R(2,401),SL(401)
      COMMON/ND/ON(401),HN(401),N2N(401),O2N(401),PHION(401),TN(401)
      COMMON/ALT/Z(401),DT,DH,THF,EPS,TF,ITF,JMAX,JMAX1,ITER,ION
      COMMON/SAV/NSAVE(2,401),TISAV(3,401),FY(2,401),UN(401),EHT(3,401)
      COMMON/APAVLOV/AKEI,OFAC,IPAVLOV
      DATA IPAVLOV/0/   !... turns on Pavlovs cooling rates if .NE.0
      DATA AKEI/1.0/    !... factor for e-i cooling (standard=1.0)
      DATA OFAC/1.0/    !... factor for i-O cooling (standard=1.0)
      IPR=3          !... signifies Temp sol
      ID=3           !... unit for writing diagnostics in Nyyyyddd.dat
      IF(J.GT.JMAX/2) ID=9    !... for Syyyyddd.dat

      !.... set up the ion and electron energy equations
      CALL SOLVE_IONHEAT(J,ILJ,ID,N,F,TI,TISAV,EHT,HE)
      CALL SOLVE_ELECTRON_HEAT(J,ILJ,ID,N,F,TI,TISAV,EHT,HE)
      RETURN
      END

C::::::::::::::::::::::::::: SOLVE_ELECTRON_HEAT :::::::::::::::::::::
C.. This subroutine sets up the integrated electron energy equation for
C.. the Newton solver F(2)
C.. It calculates the electron cooling rate, the derivative of 
C.. temperature,T, wrt to time , the heating rate, heat_flux and the
C.. loss rate for three adjacent grid points j,j-1 and j+1. These are 
C.. stored and later used to determine the values at the upper and lower
C.. midpoints by interpolation. The upper midpoint corresponds to the
C.. point midway between j and j+1 whereas lower midpoint corresponds to
C.. that between j-1 and j. The other terms are Q_UPPR and Q_LOWR which
C.. are the integrated values of the heat flow between j+1 and j && j and
C.. j-1 respectively.The basic equation is 
C.. F(2)= DT/DT(DERIVT) + DIV(HEAT_FLUX) + ELEC_ION_COOLING_RATE(E_CLRT)
C..       + ELEC_NEUT_COOLING(HLOSS_E) - HEATING_RATE(HT_RAT) = 0
C.. TERLIN interpolates and the integration is done between the midpoints.
C-- Modified in August 1995 to include ExB drift term. Bailey et al.
C-- PSS 1983, page 389 (DIVEM)

      SUBROUTINE SOLVE_ELECTRON_HEAT(JI,ILJ,ID,N,F,TI,TISAV,EHT,HE)
      IMPLICIT DOUBLE PRECISION(A-H,K-L,N-Z)
      DIMENSION N(4,401),TI(3,401),TISAV(3,401),f(20),
     > DER_T(3),CLRT(3),HLOSE(3),DEL_S(2),EHT(3,401),
     > HT_ERG(3),KE(3),DIV_VEL(3),HE(401)
      COMMON/VN/U(2,401),BG(401),BM(401),GR(2,401),R(2,401),SL(401)
      COMMON/ND/ON(401),HN(401),N2N(401),O2N(401),PHION(401),TN(401)
      COMMON/ALT/Z(401),DT,DH,FHT,EPS,TF,ITF,JMAX,JMAX1,ITER,JION
      COMMON/EXB/EQ_VPERP,DIV_VEM(401)

      DATA  BOLTZ /1.3807E-16/

      !... evaluate terms at J-1, J, J+1 for integration, JI is the
      !... actual grid point. All terms divided by magnetic field (BM)
      DO 300 I=1,3
        J=I+JI-2
        NE = N(1,J)+N(2,J)+N(3,J)+N(4,J)   !... electron density
        !... temperature derivative dT/dt 
        DER_T(I) = 1.5*BOLTZ*NE*(TI(3,J)-TISAV(3,J))/BM(J)/DT
        !... electron heating rate in ergs
        HT_ERG(I) = EHT(3,J) * 1.6E-12/BM(J)
        !.. determine the thermal conductivity KE 
        CALL GET_THERMCONDUCTIVITY(TI,NE,HE,KE,I,J)
        !.. ion and neutral cooling rates
        CALL GET_LOSS(I,J,TI,N,BM,HLOSE,CLRT)
        !.. ExB drift term from Bailey's eqns. 16,19 in PSS,1983,p389
        DIV_VEL(I) = BOLTZ*DIV_VEM(J) *NE * TI(3,J)/BM(J)
 300  CONTINUE

      !.. distances along field line between grid points for the integration 
      DEL_S(2) = SL(JI+1) - SL(JI)
      DEL_S(1) = SL(JI) - SL(JI-1)

      !... Now perform the integration. DERIVT, E_CLRT, LOSS_E, and 
      !... HT_RAT are the four values after integration 
      CALL TERLIN(JI,DER_T(1),DER_T(2),DER_T(3),DEL_S(1),DEL_S(2),
     >  DTL,DERIVT,DTU)
      CALL TERLIN(JI,CLRT(1),CLRT(2),CLRT(3),DEL_S(1),DEL_S(2),
     >  CLRT_L,E_CLRT,CLRT_U)
      CALL TERLIN(JI,HLOSE(1),HLOSE(2),HLOSE(3),DEL_S(1),DEL_S(2),
     >  HLOSE_L,LOSS_E,HLOSE_U)
      CALL TERLIN(JI,HT_ERG(1),HT_ERG(2),HT_ERG(3),DEL_S(1),DEL_S(2),
     >  EH_L,HT_RAT,EH_U)
      CALL TERLIN(JI,DIV_VEL(1),DIV_VEL(2),DIV_VEL(3),DEL_S(1),DEL_S(2),
     >  DRFT_L,DRFT_TRM,DRFT_U)

      !.... Calculate the divergence of the heat flux term
      KU_AVG= (KE(3) + KE(2))/2.0
      KL_AVG = (KE(2) + KE(1))/2.0
      BU_AVG = (BM(JI) + BM(JI+1)) * 0.5
      BL_AVG = (BM(JI) + BM(JI-1)) * 0.5
      Q_UPPR = KU_AVG *(TI(3,JI+1)-TI(3,JI))/DEL_S(2)/BU_AVG
      Q_LOWR = KL_AVG *(TI(3,JI)-TI(3,JI-1))/DEL_S(1)/BL_AVG

      !... now form the error function for Te
      F(2) = DERIVT + E_CLRT + LOSS_E - HT_RAT -  Q_UPPR + Q_LOWR
     >  + DRFT_TRM*EQ_VPERP
 
      !... This section for diagnostic print changed 2022-10-15
      IF(ILJ.EQ.4) THEN
        WRITE(ID,501)
        WRITE(ID,5555) Z(JI),NINT(TI(3,J)),NINT(TI(2,J)),NINT(TN(J))
     >   ,F(2),DERIVT,E_CLRT,LOSS_E,Q_UPPR,Q_LOWR
     >   ,HT_RAT,DRFT_TRM*EQ_VPERP,Q_UPPR-Q_LOWR,NE
      ENDIF
 501  FORMAT('     ALT    Te     Ti    Tn      F       dT/dt'
     >   ,7X,'C_ei       C_en      Q_UPPR     Q_LOWR      H_tot'
     >   ,5X,'DIV_VP      DivQ        Ne')
 5555 FORMAT(F9.1,3I6,1P,22E11.3)
      RETURN
      END
 
C::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
C   ROUTINE:	GET_LOSS
C   PURPOSE:	calculates the cooling rates, hlose for electrons
C   ARGUMENT LIST: 
C   REFERENCES: Rees and Roble, Rev. Geophys.(1975) P220 
C               Schunk and Nagy Rev Geophys (1978) p366
C   NAME		USE
C   ----		---		
C   I,J			LOOP VARIABLES
C   TI			TEMPERATURE
C   NE			ELECTRON DENSITY
C   HLOSE		loss OF ELECTRONS
C   BM			MAGNETIC FIELD  			
C::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
      SUBROUTINE GET_LOSS(I,J,TI,N,BM,HLOSE,CLRT)
      IMPLICIT DOUBLE PRECISION(A-H,K-L, N-Z)
      REAL PAVLOV_VN2,PAVLOV_VO2,PAVLOV_OFS
      DIMENSION HLOSE(3), TI(3,401), BM(401),CLRT(3),N(4,401)
      COMMON/ND/ON(401),HN(401),N2N(401),O2N(401),PHION(401),TN(401)
      COMMON/ALT/Z(401),DT,DH,FHT,EPS,TF,ITF,JMAX,JMAX1,ITER,JION
      COMMON/APAVLOV/AKEI,OFAC,IPAVLOV

      !.... electron-ion cooling rate
      NE = N(1,J)+N(2,J)+N(3,J)+N(4,J)
      NI = N(1,J)/16 + N(2,J) + N(3,J)/30.0

      !.. Electron- ion cooling rate KEI. The Coulomb logarithm from Itikawa
      !..  JATP 1975, page 1001. Previously assumed to be 15.0
      DEBYE = 1.602E-12 * SQRT(4.0*3.1415926*NE/(1.381E-16*TI(3,J)))
      COULOG = 0.43086 + DLOG(TI(3,J)) - DLOG(DEBYE)
      CLRT(I) = 1.232e-17*NE*NI*(TI(3,J)-TI(2,J))/BM(J)/(TI(3,J)**1.5)
      CLRT(I)= CLRT(I)*COULOG/15.0               !.. Itikawa correction
      CLRT(I)= AKEI*CLRT(I)                      !.. fudge factor

      !... Neglect losses to neutrals above 1000 km
      IF(O2N(J).LT.1.0E-10) THEN
         HLOSE(I)=0.0
         RETURN
      ENDIF

      !.. These terms are evaluated here for efficiency because
      !.. they are used in many places
      TE2=TI(3,J)
      RTE2=1.0/TE2
      TNJ=TN(J)
      RTNJ=1.0/TNJ
      SQTE=SQRT(TE2)
      TDIF=TE2-TNJ
      TDIFRT=TDIF*RTE2*RTNJ

       !..  elec-neut elas coll loss rates s&N(1978) rev geophys p366
       LEN2=1.77E-19*NE*N2N(J)*(1.0-1.2E-4*TE2)*TE2*TDIF
       LEO2=1.21E-18*NE*O2N(J)*(1.0+3.6E-2*SQTE)*SQTE*TDIF
       LEO=7.9E-19*NE*ON(J)*(1.0+5.7E-4*TE2)*SQTE*TDIF
       LEH=9.63E-16*NE*HN(J)*(1.0-1.35E-4*TE2)*SQTE*TDIF

       !...  rotational loss rates b&k 268 ####
       LRN2=2.0E-14*NE*N2N(J)*TDIF/SQTE
       LRO2=7.0E-14*NE*O2N(J)*TDIF/SQTE

       !...  n2 vib loss rates b&k p268 and s&n p364  +++++
       EF=1.06E+4+7.51E+3*TANH(1.1E-3*(TE2-1800.0))
       GE=3300.0+1.233*(TE2-1000.0)-2.056E-4*(TE2-1000.0)*
     >   (TE2-4000.0)
       LVN2=-2.99E-12*NE*N2N(J)*EXP(5.0E-4*RTE2*EF*(TE2-2000.0))
     > *(EXP(-GE*TDIFRT)-1.0)

       !... o2 vib loss rate s&n p364  $$$$
       HS=3300.0-839.0*SIN(1.91E-4*(TE2-2700.0))
       LVO2=-5.196E-13*NE*O2N(J)*EXP(1.4286E-3*RTE2*HS*(TE2-700.0))
     >        *(EXP(-2770.0*TDIFRT)-1.0)

       !... fine structure cooling by atomic oxygen
       D1=EXP(-228.0*RTNJ)
       D2=EXP(-326.0*RTNJ)
       E1=EXP(-228.0*RTE2)
       E2=EXP(-326.0*RTE2)
       E3=EXP(-98.0*RTE2)*D1
       LF1=8.49E-6*TE2**0.519*(0.02*(D1-E1)-5.91E-9*
     >     TDIF*(2.019*D1+(228.0*RTE2+2.019)*E1))
       LF2=7.7E-6*TE2**0.3998*(0.028*(D2-E2)-5.91E-9*
     >     TDIF*(1.8998*D2+(326.0*RTE2+1.8998)*E2))
       LF3=2.22E-7*TE2**0.768*(0.008*(D2-E3)-5.91E-9*
     >     TDIF*(2.268*D2+(98.0*RTE2+2.268)*E3))
       ZFO=5.0+3.0*D1+D2
       LFO=-8.629E-6*NE*ON(J)*(LF1+LF2+LF3)/ZFO
 
       !.. excitation of O(1D)  S&N P365. 
       DE=2.4E4+(TE2-1500.0)*(0.3-1.947E-5*(TE2-4000.0))
       EXH=3.3333E-4*DE*(TE2-3000.0)*RTE2
       IF(EXH.GT.70.0)EXH=70.0
       EXH2=22713.0*TDIFRT
       IF(EXH2.GT.70.0)EXH2=70.0
       LF1D=-1.57E-12*NE*ON(J)*EXP(EXH)*(EXP(-EXH2)-1.0)


         !.. Electron heating from N2v
         !.. RLVN2=1.0
         !.. IF(DABS(LVN2).GT..0001) RLVN2=PLVN2/LVN2
         !.. LVN2=PLVN2
         !.. IF(RLVN2.GT.2.OR.RLVN2.LT.0.5) 
         !.. WRITE(6,'(2F6.0,1P,9E11.2)') TE2,TNJ,PLVN2,LVN2
         !... lvn2=1.0*lvn2*RVBT(J)  !... for N2*

       !.. N2*, O2*, Ofs from Pavlov. 
       IF(IPAVLOV.NE.0) THEN
          IF(IPAVLOV.GT.0) WRITE(6,*) ' PAVLOV electron cooling'
          IF(IPAVLOV.GT.0) IPAVLOV=-99
          LVN2=PAVLOV_VN2(TE2,TNJ)*NE*N2N(J)
          LVO2=PAVLOV_VO2(TE2,TNJ)*NE*O2N(J)
          LFO =PAVLOV_OFS(TE2,TNJ)*NE*ON(J)
       ENDIF

 14    KEN=LEN2+LEO2+LEO+LEH+LRN2+LRO2
       HLOSE(I)=(KEN+LVN2+LVO2+LFO+LF1D)*1.6E-12/BM(J)
       RETURN
       END

C::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
C     ROUTINE:	GET_THERMCONDUCTIVITY
C.....PURPOSE:	calculates ke, the thermal conductivity co-efficient
C
C::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: 
      SUBROUTINE GET_THERMCONDUCTIVITY(TI,NET,HE,KE,I,J)
      IMPLICIT DOUBLE PRECISION(A-H,K-L,N-Z)
      DIMENSION KE(3),TE(3), TI(3,401),HE(401)
      COMMON/ND/ON(401),HN(401),N2N(401),O2N(401),PHION(401),TN(401)

C.....TE is the electron temperature
      TE(I)=TI(3,J)
      SQTE=SQRT(TE(I))
      KN2=(2.82E-17*SQTE-3.41E-21*SQTE*TE(I))*N2N(J)
      KO2=(2.2E-16+7.92E-18*SQTE)*O2N(J)
      KO=3.4E-16*ON(J)
      KHE=5.6E-16*HE(J)
      KH=(5.47E-15-7.45E-19*TE(I))*HN(J)
      !.. KN is Bank's factor for neutral effects on conductivity.
      !.. 1.0E+4 is added to avoid numerical problems for low e densities
      KN=3.22E+4*TE(I)**2*(KN2+KO2+KO+KHE+KH)/(NET + 1.0E+4)
      KE(I)=1.232E-6*(TE(I)**2.5)/(1.0+KN)   !-- Thermal conductivity 
      RETURN
      END

C:::::::::::::::::::::::::: SOLVE_IONHEAT ::::::::::::::::::::::::::
C.. This subroutine is similar to SOLVE_IONHEAT. Look at the Te function
C.. for more explanation. It determines the interpolated production and
C.. loss processes for ions. It calls RATES to get the rate constants
C.. and TERLIN to do the interpolation
      SUBROUTINE SOLVE_IONHEAT(JI,ILJ,ID,N,F,TI,TISAV,EHT,HE)
      IMPLICIT DOUBLE PRECISION(A-H,K-L,N-Z)
      DIMENSION N(4,401),TI(3,401),TISAV(3,401),F(20),EHT(3,401)
     > ,DER_T(3),CLRT(3),HLOSI(3),KI(3), DEL_S(2),DIV_VEL(3)
     > ,HT_ERG(3),HE(401),COOL(22)
      COMMON/VN/U(2,401),BG(401),BM(401),GR(2,401),R(2,401),SL(401)
      COMMON/ND/ON(401),HN(401),N2N(401),O2N(401),PHION(401),TN(401)
      COMMON/ALT/Z(401),DT,DH,FHT,EPS,TF,ITF,JMAX,JMAX1,ITER,JION
      COMMON/DN2/TEJ(401),TIJ(401),NUX(2,401),RBM(401),UNJ(401),DS(401)
      COMMON/EXB/EQ_VPERP,DIV_VEM(401)
      COMMON/APAVLOV/AKEI,OFAC,IPAVLOV

      DATA BOLTZ/1.3807E-16 /
  
      !... evaluate terms at j-1, j, J+1 for integration, JI is the
      !... actual grid point. All terms divided by magnetic field (BM)
      DO 300 I=1,3
        J=I+JI-2
        TNJ = TN(J)
        NE = N(1,J) + N(2,J) + N(3,J)+N(4,J)
        NI = N(1,J)/16 + N(2,J) + N(3,J)/30.0
        RBM(J) = 1/BM(J)
        DER_T(I) = 1.5*BOLTZ*NE*(TI(2,J) - TISAV(2,J))*RBM(J)/DT
        !.. direct heating rate for ions (ring current?)
        HT_ERG(I) = EHT(1,J) * 1.6E-12/BM(J)

        !.. electron-ion heating/cooling rate
        !.. Electron- ion cooling rate KEI. The Coulomb logarithm from Itikawa
        !..  JATP 1975, page 1001. Previously assumed to be 15.0
        DEBYE = 1.602E-12 * SQRT(4.0*3.1415926*NE/(1.381E-16*TI(3,J)))
        COULOG = 0.43086 + DLOG(TI(3,J)) - DLOG(DEBYE)
        CLRT(I)=1.232E-17*NE*NI*(TI(3,J)-TI(2,J))*RBM(J)/(TI(3,J)**1.5)
        CLRT(I)= CLRT(I)*COULOG/15.0               !.. Itikawa correction
        CLRT(I)= AKEI*CLRT(I)                      !.. fudge factor
    

        !.. loss rate coeff to neutrals (kin) rees & roble(1975) P220 ,,,
        !.. RCE = resonant charge exchange, pol = polarization interaction
        RCE=(0.21*N(1,J)*ON(J)*OFAC+1.4*N(2,J)*HN(J))
        POL1=N(1,J)*(6.6*N2N(J)+2.8*HE(J)+5.8*O2N(J)+5.6*ON(J)*OFAC)
     >       +N(2,J)*(5.5*HE(J)+1.9*HN(J))
        POL2=(0.36*N(1,J)*HN(J)+0.4*N(2,J)*ON(J))*SQRT(TNJ)
        KIN=(RCE*SQRT(TNJ+TI(2,J))+POL1+POL2)*1.6E-26
        CALL ION_NEUTRAL(0,N(1,J),N(2,J),ON(J),N2N(J),O2N(J)
     >   ,HN(J),HE(J),TN(J),TI(2,J),Z(J),COOL)

        HLOSI(I)=+KIN*(TI(2,J)-TNJ)*RBM(J)
        HLOSI(I)=+COOL(1)*(TI(2,J)-TNJ)*RBM(J)

        !.. determine the thermal conductivity KI 
        KI(I)=1.84E-8*(N(1,J)+4*N(2,J))*(TI(2,J)**2.5)/NE
        DIV_VEL(I) = DIV_VEM(J) * BOLTZ* NE*TI(2,J)/BM(J)
 300  CONTINUE

       !.. distances between grid points
       DEL_S(2) = SL(JI+1) - SL(JI)
       DEL_S(1) = SL(JI) - SL(JI-1)

      !.. Do interpolations. derivT, E_CLRT, loss , drft_trm are the
      !.. values after integration 
      CALL TERLIN(JI,DER_T(1),DER_T(2),DER_T(3),DEL_S(1),DEL_S(2),
     >	DTL,DERIV_T,DTU)
      CALL TERLIN(JI,CLRT(1),CLRT(2),CLRT(3),DEL_S(1),DEL_S(2),
     >  CLRT_L,E_CLRT,CLRT_U)
      CALL TERLIN(JI,HLOSI(1),HLOSI(2),HLOSI(3),DEL_S(1),DEL_S(2),
     >  HLOSI_L,LOSS,HLOSI_U)
      CALL TERLIN(JI,DIV_VEL(1),DIV_VEL(2),DIV_VEL(3),DEL_S(1),DEL_S(2),
     >  DRFT_L,DRFT_TRM,DRFT_U)
      CALL TERLIN(JI,HT_ERG(1),HT_ERG(2),HT_ERG(3),DEL_S(1),DEL_S(2),
     >  EH_L,HT_RAT,EH_U)

      !.. ion thermal conductivity
      KU_AVG= (KI(3) + KI(2))/2.0
      KL_AVG = (KI(2) + KI(1))/2.0
      TU_AVG = (TI(2,JI+1)+TI(2,JI))/2.0
      TL_AVG = (TI(2,JI)+TI(2,JI-1))/2.0
      BU_AVG = (BM(JI) + BM(JI+1)) * 0.5
      BL_AVG = (BM(JI) + BM(JI-1)) * 0.5
      Q_UPPR = KU_AVG *(TI(2,JI+1)-TI(2,JI))/DEL_S(2)/BU_AVG
      Q_LOWR = KL_AVG *(TI(2,JI)-TI(2,JI-1))/DEL_S(1)/BL_AVG

      !.. integrated energy equation for ions
      F(1) = DERIV_T-E_CLRT+LOSS-Q_UPPR+Q_LOWR+DRFT_TRM*EQ_VPERP-HT_RAT
      
      !.. diagnostic print
      IF(ILJ.EQ.5) THEN
        WRITE(ID,501)
        WRITE(ID,5555) Z(JI),F(1),DERIV_T,E_CLRT,LOSS,Q_UPPR,Q_LOWR
     >      ,DRFT_TRM*EQ_VPERP
      ENDIF
 501  FORMAT(' Ti: ALT     F        dT/dt      C_ie      C_in'
     >   ,6X,'Q_UPPR   Q_LOWR    DIV_VP')
 5555   FORMAT(F9.1,1P,22E10.2)
      RETURN
      END
C::::::::::::::::::::::: ECRATS :::::::::::::::::::::::::::::::::::
C....... Subroutine for printing thermal electron cooling rates
C....... This function needs to be eliminated eventually and replaced
C....... by GET_loss. It is used to print cooling rates in the Hyyyyddd
C....... file of the FLIP run and also in the FLIPRINT phase 
      SUBROUTINE ECRATS(ID,J,ALT,TE,TI,TN,OPLS,HPLS,NMIN,OX,N2,O2,HN
     > ,HE,EHT,TLSS,LFO,COOL)
      IMPLICIT DOUBLE PRECISION(A-H,K,L,N,O-Z)
      REAL PAVLOV_VN2,PAVLOV_VO2,PAVLOV_OFS
      COMMON/V91/RVBT(401)
      COMMON/APAVLOV/AKEI,OFAC,IPAVLOV
      DIMENSION FD(9),COOL(22),COOLI(22) !.. heating and cooling rates
C
      IDID=ID
      IF(IDID.EQ.7)IDID=9
      IF(ID.NE.IDS.AND.ID.NE.0) WRITE(IDID,757)
      IDS=ID
 757  FORMAT(/,8X,'Thermal electron cooling(C_) and heating(H_) rates'
     > /,6X,'ALT  Te   Ti   Tn     C_ei      C_rN2     C_rO2    C_vN2',
     > '    C_vO2     C_fsO      C_O1D     C_tot      H_tot     H_N2D'
     > '     HOp2D     TIC       ROPOX     RHPHN     POPN2     POPHE'
     > '     POPO2     POPON     PHPHE     PHPHN     POPHN     PHPON')
        TLSS=0.0
        LFO=0.0

        !.. Don't calculate cooling rates if neutral densities too low
        IF(O2.LT.1.0E-10) RETURN     
        SQTE=SQRT(TE)
        NE=OPLS+HPLS+NMIN
        NI=OPLS/16.0+HPLS+NMIN/30.0

        !.. Electron- ion cooling rate KEI. The Coulomb logarithm from Itikawa
        !..  JATP 1975, page 1001. Previously assumed to be 15.0
        DEBYE = 1.602E-12 * SQRT(4.0*3.1415926*NE/(1.381E-16*TE))
        COULOG = 0.43086 + DLOG(TE) - DLOG(DEBYE)
        KEI=1.232E-17*NE*NI*(TE-TI)/(TE**1.5)/1.6E-12    !.. e - i cooling
        KEI= KEI*COULOG/15.0                     !.. Itikawa correction
        KEI= AKEI*KEI                            !.. fudge factor


        RTE2=1.0/TE
        RTNJ=1.0/TN
 13     TDIF=TE-TN
        TDIFRT=TDIF*RTE2*RTNJ
        LEN2=1.77E-19*NE*N2*(1.0-1.2E-4*TE)*TE*TDIF
        LEO2=1.21E-18*NE*O2*(1.0+3.6E-2*SQTE)*SQTE*TDIF
        LEO=7.9E-19*NE*OX*(1.0+5.7E-4*TE)*SQTE*TDIF
        LEH=9.63E-16*NE*HN*(1.0-1.35E-4*TE)*SQTE*TDIF
C
C###  ROTATIONAL loss RATES B&K 268 ####
        LRN2=2.0E-14*NE*N2*TDIF/SQTE
        LRO2=7.0E-14*NE*O2*TDIF/SQTE
C
C+++  N2 VIB loss RATES B&K P268 AND S&N P364  +++++
        EF=1.06E+4+7.51E+3*TANH(1.1E-3*(TE-1800.0))
        GE=3300.0+1.233*(TE-1000.0)-2.056E-4*(TE-1000.0)*(TE-4000.0)
        EXH=5.0E-4*RTE2*EF*(TE-2000.0)
        IF(EXH.GT.70.0)EXH=70.0
        EXH2=GE*TDIFRT
        IF(EXH2.GT.70.0)EXH2=70.0
        LVN2=-2.99E-12*NE*N2*EXP(EXH)
     > *(EXP(-EXH2)-1.0)
C...... LVN2=1.0*LVN2*RVBT(J)
C
C$$$$ O2 VIB loss RATE S&N P364  $$$$
        HS=3300.0-839.0*SIN(1.91E-4*(TE-2700.0))
        EXH=1.4286E-3*RTE2*HS*(TE-700.0)
        IF(EXH.GT.70.0)EXH=70.0
        EXH2=2770.0*TDIFRT
        IF(EXH2.GT.70.0)EXH2=70.0
        LVO2=-5.196E-13*NE*O2*EXP(EXH2)
     >        *(EXP(-EXH2)-1.0)
C
C;;;;;  FINE STRUCTURE EXCITATIONS S&N P365 ;;;;;;
C;;;;; NOTE THAT D1,D2,E1 ETC. MAY NEED TO BE CHANGED
C;;;; NOTE THAT A TERM IS ADDED TO NE AT LOW ALTS FOR O2+ , NO+
        HD1=228.0*RTNJ
        IF(HD1.GT.70.0)HD1=70.0
        D1=EXP(-HD1)
        HD2=326.0*RTNJ
        IF(HD2.GT.70.0)HD2=70.0
        D2=EXP(-HD2)
        HE1=228.0*RTE2
        IF(HE1.GT.70.0)HE1=70.0
        E1=EXP(-HE1)
        HE2=326.0*RTE2
        IF(HE2.GT.70.0)HE2=70.0
        E2=EXP(-HE2)
        HE3=98.0*RTE2
        IF(HE3.GT.70.0)HE3=70.0
        E3=EXP(-HE3)*D1
        LF1=8.49E-6*TE**0.519*(0.02*(D1-E1)-5.91E-9*
     >     TDIF*(2.019*D1+(228.0*RTE2+2.019)*E1))
        LF2=7.7E-6*TE**0.3998*(0.028*(D2-E2)-5.91E-9*
     >     TDIF*(1.8998*D2+(326.0*RTE2+1.8998)*E2))
        LF3=2.22E-7*TE**0.768*(0.008*(D2-E3)-5.91E-9*
     >     TDIF*(2.268*D2+(98.0*RTE2+2.268)*E3))
        ZFO=5.0+3.0*D1+D2
        LFO=-8.629E-6*NE*OX*(LF1+LF2+LF3)/ZFO
C....... LFO=0.333*LFO
C++++  EXCITATION OF O(1D) BY PE'S S&N P365. - SMALL ONLY CALC FOR PRINT
        DE=2.4E4+(TE-1500.0)*(0.3-1.947E-5*(TE-4000.0))
        EXH=3.3333E-4*DE*(TE-3000.0)*RTE2
        IF(EXH.GT.70.0)EXH=70.0
        EXH2=22713.0*TDIFRT
        IF(EXH2.GT.70.0)EXH2=70.0
        LF1D=-1.57E-12*NE*OX*EXP(EXH)*(EXP(-EXH2)-1.0)

       IF(IPAVLOV.NE.0) THEN
          LVN2=PAVLOV_VN2(TE,TN)*NE*N2
          LVO2=PAVLOV_VO2(TE,TN)*NE*O2
          LFO =PAVLOV_OFS(TE,TN)*NE*OX
       ENDIF

C......... add other losses to electron-ion
      TLSS=LEN2+LEO2+LEO+LEH+LRO2+LRN2+LVO2+LVN2+LF1D
      IF(TLSS.LT.0) TLSS=0.0

C......... electron heating from N(2D) .....
      DO 1199 III=1,8
 1199 FD(III)=0.0
      !..IF(ID.NE.0.AND.ID.NE.3) 
      CALL CMINOR(0,J,0,IPR,FD,3,HE)
      CALL ION_NEUTRAL(0,OPLS,HPLS,OX,N2,O2,HN,HE,TN,TI,ALT,COOLI)
      IF(ID.NE.0) WRITE(IDID,54) ALT,INT(TE),INT(TI),INT(TN),KEI
     > ,LRN2,LRO2,LVN2,LVO2,LFO,LF1D,TLSS+LFO+KEI,EHT,FD(8)*2.4,
     > FD(6)*3.31,(COOLI(III)*(TI-TN)/1.6E-12,III=1,11)
 54    FORMAT(F9.1,3I5,1P,33E10.2)

      COOL(1)=EHT
      COOL(2)=KEI
      COOL(3)=LVN2
      COOL(4)=LFO+LF1D
      COOL(5)=TLSS+LFO+KEI
      RETURN
      END
C:::::::::::::::::::: PAVLOV_VN2 :::::::::::::::::::::::::::::::::::
C... This function returns the coefficient Q0v for electron cooling by
C... vibrational excitation of N2
C... Q0v is used in equation 11. Q1v is ignored because small
C... See Pavlov Ann. Geophysicae, 1998 page 176
      REAL FUNCTION PAVLOV_VN2(TE,TN)
      IMPLICIT NONE
      INTEGER I
      REAL SUM,E1,A0LTE,B0LTE,C0LTE,D0LTE,F0LTE,AQLTE
      REAL A0(10),B0(10),C0(10),D0(10),F0(10),ALQ0(10)
      DOUBLE PRECISION TE,TN
      !... coefficients for Te < 1500
      DATA A0LTE/-6.462/,B0LTE/3.151E-2/,C0LTE/-4.075E-5/
      DATA D0LTE/2.439E-8/,F0LTE/-5.479E-12/, E1/3353./

      !... coefficients for Te >= 1500
      DATA A0/2.025,-7.066,-8.211,-9.713,-10.353,-10.819,-10.183,
     >  -12.698,-14.710,-17.538/
      DATA B0/8.782E-4,1.001E-2,1.092E-2,1.204E-2,1.243E-2,1.244E-2,
     >  1.185E-2,1.309E-2,1.409E-2,1.600E-2/
      DATA C0/2.954E-7,-3.369E-6,-3.369E-6,-3.732E-6,-3.850E-6,
     >  -3.771E-6,-3.570E-6,-3.952E-6,-4.249E-6,-4.916E-6/
      DATA D0/-9.562E-11,4.436E-10,4.891E-10,5.431E-10,5.600E-10,
     >  5.385E-10,5.086E-10,5.636E-10,6.058E-10,7.128E-10/
      DATA F0/7.252E-15,-2.449E-14,-2.706E-14,-3.008E-14,-3.1E-14,
     >  -2.936E-14,-2.769E-14,-3.071E-14,-3.3E-14,-3.941E-14/

      SUM=0.0
      !... Low temperature expression
      IF(TE.LT.1500) THEN
         AQLTE=A0LTE+B0LTE*TE+C0LTE*TE**2.+D0LTE*TE**3.+
     >       F0LTE*TE**4.-16
         SUM=(10**AQLTE)*(1-EXP(E1*(1.0/TE-1.0/TN)))

      ELSE
         !... For Te >= 1500
         DO I=1,10
           ALQ0(I)=A0(I)+B0(I)*TE+C0(I)*TE**2.+D0(I)*TE**3.+
     >         F0(I)*TE**4.-16.
           SUM = SUM+(10**ALQ0(I))*(1.-EXP(REAL(I)*E1*(1.0/TE-1.0/TN)))
         ENDDO

      ENDIF

      !... return the coefficient Q0v
      PAVLOV_VN2=(1-EXP(-E1/TN))*SUM
      RETURN
      END
C::::::::::::::::: PAVLOV_VO2 :::::::::::::::::::::::::::::::::
C.... This function is Pavlov's electron cooling rate by
C.... vibrational excitation of O2. 
C.... Reference: Ann. Geophys. 16, 1007, 1998
      REAL FUNCTION PAVLOV_VO2(TE,TN)
      IMPLICIT NONE
      INTEGER I
      REAL E1,A(7),B(7),C(7),D(7),F(7),G(7),Q(7),SUM
      DOUBLE PRECISION TE,TN
      DATA A/8.56E-15,1.15E-18,1.77E-23,7.05E-25,2.14E-28,2.94E-31,
     > 1.08E-35/
      DATA B/303,305,301,301,299,300,299/
      DATA C/10.,18.78,29.,30.07,37.52,43.27,53.1/
      DATA D/-0.2,-0.25,-0.2,-0.31,-0.42,-0.39,-0.9/
      DATA F/1.05E-3,9.24E-4,6.16E-4,1.00E-3,5.28E-4,7.85E-4,5.7E-4/
      DATA G/3150,3450,3150,1800,3000,3000,2900/
      DATA E1/2239./

      SUM=0.0
      DO I=1,7
        Q(I)=A(I)*EXP((1-B(I)/TE)*(C(I)+D(I)*SIN(F(I)*(TE-G(I)))))
        SUM = SUM+Q(I)*(1.-EXP(I*E1*(1./TE-1./TN)))
      ENDDO

      PAVLOV_VO2=SUM
      RETURN
      END
C:::::::::::::::::: PAVLOV_OFS :::::::::::::::::::::::::::
C.... This is the O fine structure cooling rate of Pavlov
C.... Ann. Geophys. 17, 919, 1999
      REAL FUNCTION PAVLOV_OFS(TE,TN)
      IMPLICIT NONE
      REAL S21,S20,S10,D,FTN
      DOUBLE PRECISION TE,TN
      S21=1.863E-11
      S20=1.191E-11
      S10=8.249E-16*TE**0.6*EXP(-227.7/TN)
      D=5.0+EXP(-326.6/TN)+3*EXP(-227.7/TN)
      FTN=8.132E-9/D
      PAVLOV_OFS=FTN*(TE-TN)/(TE*TN)
      RETURN
      END
C:::::::::::::::: ION_NEUTRAL ::::::::::::::::::::::
C.... loss rate coeff to neutrals (kin) rees & roble(1975) P220 ,,,
      SUBROUTINE ION_NEUTRAL(IWR,OPLUS,HPLUS,ON,N2N,O2N,HN,HEN
     >   ,TN,TI,ALT,COOL)
      IMPLICIT NONE
      INTEGER IWR,IPAVLOV,IK
      DOUBLE PRECISION OPLUS,HPLUS,ON,N2N,O2N,HN,HEN,TN,TI,KIN,OFAC
     > ,RCE,POL1,POL2,ALT2,AKEI,ALT,COOL(22)
      COMMON/APAVLOV/AKEI,OFAC,IPAVLOV
        !.. Resonant charge exchange
        COOL(2)=1.6E-26*(0.21*OPLUS*ON*OFAC)*SQRT(TI+TN)
        COOL(3)=1.6E-26*(1.4*HPLUS*HN)*SQRT(TI+TN)
        !... Polarization interactions
        COOL(4)=1.6E-26*OPLUS*(6.6*N2N)
        COOL(5)=1.6E-26*OPLUS*(2.8*HEN)
        COOL(6)=1.6E-26*OPLUS*(5.8*O2N)
        COOL(7)=1.6E-26*OPLUS*(5.6*ON*OFAC)
        COOL(8)=1.6E-26*HPLUS*(5.5*HEN)
        COOL(9)=1.6E-26*HPLUS*(1.9*HN)
        COOL(10)=1.6E-26*(0.36*OPLUS*HN)*SQRT(TN)
        COOL(11)=1.6E-26*(0.40*HPLUS*ON)*SQRT(TN)
        !... total cooling rate
        COOL(1)=COOL(2)+COOL(3)+COOL(4)+COOL(5)+COOL(6)+COOL(7)+
     >    COOL(8)+COOL(9)+COOL(10)+COOL(11)
      RETURN
      END
